#include "AnaDef.h"
#include "style.h"
#include "nevent.h"

TString fTag;
Double_t PsN = -999;
Double_t PsVrecN = -999;
Double_t PsVrecVzN = -999;

void ini(const TString fin)
{
  Int_t kcen = -999;
  if(fin.Contains("CEN0"))
    kcen = 0;
  else if(fin.Contains("CEN1"))
    kcen = 1;
  else if(fin.Contains("CEN2"))
    kcen = 2;
  else if(fin.Contains("CEN3"))
    kcen = 3;
  else if(fin.Contains("CEN4"))
    kcen = 4;
  else if(fin.Contains("CEN5"))
    kcen = 5;
  else if(fin.Contains("CEN6"))
    kcen = 6;

  //===

  fTag = fin;

  fTag.ReplaceAll("list.","");
  fTag=fTag(0, fTag.First('.'));
  if(fTag.Contains("_"))
    fTag=fTag(0, fTag.First('_'));

  fTag=fTag(fTag.First('-')+1,fTag.Length());
  fTag=fTag(0,fTag.First('-'));

  if(kcen!=-999)
    fTag+=Form("CEN%d", kcen);

  if(!fTag.Contains("MC"))
    fTag+="DATA";

  //===========
  if(fTag.Contains("10dMC")){
    PsN = PsNLHC10dMC;
    PsVrecN = PsVrecNLHC10dMC;
    PsVrecVzN = PsVrecVzNLHC10dMC;
  }
  else if(fTag.Contains("10eMC")){
    PsN = PsNLHC10eMC;
    PsVrecN = PsVrecNLHC10eMC;
    PsVrecVzN = PsVrecVzNLHC10eMC;
  }
  else if(fTag.Contains("11aMC")){
    PsN = PsNLHC11aMC;
    PsVrecN = PsVrecNLHC11aMC;
    PsVrecVzN = PsVrecVzNLHC11aMC;
  }
  else if(fTag.Contains("10hMCCEN0")){
    PsN = PsNLHC10hMCCEN0;
    PsVrecN = PsVrecNLHC10hMCCEN0;
    PsVrecVzN = PsVrecVzNLHC10hMCCEN0;
  }
  else if(fTag.Contains("10hMCCEN1")){
    PsN = PsNLHC10hMCCEN1;
    PsVrecN = PsVrecNLHC10hMCCEN1;
    PsVrecVzN = PsVrecVzNLHC10hMCCEN1;
  }
  else if(fTag.Contains("10hMCCEN2")){
    PsN = PsNLHC10hMCCEN2;
    PsVrecN = PsVrecNLHC10hMCCEN2;
    PsVrecVzN = PsVrecVzNLHC10hMCCEN2;
  }
  else if(fTag.Contains("10hMCCEN3")){
    PsN = PsNLHC10hMCCEN3;
    PsVrecN = PsVrecNLHC10hMCCEN3;
    PsVrecVzN = PsVrecVzNLHC10hMCCEN3;
  }
  else if(fTag.Contains("10hMCCEN4")){
    PsN = PsNLHC10hMCCEN4;
    PsVrecN = PsVrecNLHC10hMCCEN4;
    PsVrecVzN = PsVrecVzNLHC10hMCCEN4;
  }
  else if(fTag.Contains("10hMCCEN5")){
    PsN = PsNLHC10hMCCEN5;
    PsVrecN = PsVrecNLHC10hMCCEN5;
    PsVrecVzN = PsVrecVzNLHC10hMCCEN5;
  }
  else if(fTag.Contains("10hMCCEN6")){
    PsN = PsNLHC10hMCCEN6;
    PsVrecN = PsVrecNLHC10hMCCEN6;
    PsVrecVzN = PsVrecVzNLHC10hMCCEN6;
  }
  else {
    printf("wrong fTag %s\n", fTag.Data()); exit(1);
  }
}

void normal(TH1D *hh, Double_t nor=-999)
{
  if(nor<0){
    nor = hh->Integral(0,hh->GetNbinsX()+1);
    if(nor<1)
      return;
  }

  for(Int_t ii=0; ii<=hh->GetNbinsX()+1; ii++){
    const Double_t cc = hh->GetBinContent(ii);
    const Double_t ee = hh->GetBinError(ii);
    const Double_t ww = hh->GetBinWidth(ii);
    hh->SetBinContent(ii, cc/nor/ww);
    hh->SetBinError(ii, ee/nor/ww);
  }
}


//up is smaller than down
void ProperDivide(TH1D *hup, const TH1D *hdown)
{
  //in normal both content and error are divided by a constant number nor*binwidth
  for(Int_t ii=0; ii<=hup->GetNbinsX()+1; ii++){
    //up
    const Double_t yup = hup->GetBinContent(ii);
    if(yup<EPSILON)
      continue;

    //y = n*scale
    //e = sqrt(n)*scale
    const Double_t eup = hup->GetBinError(ii);
    const Double_t scaleup = eup*eup/yup;
    const Double_t nup = yup/scaleup;

    const Double_t ydown = hdown->GetBinContent(ii);
    if(ydown<EPSILON){
      printf("ydown < eps! %e\n", ydown);exit(1);
    }
    const Double_t edown = hdown->GetBinError(ii);
    const Double_t scaledown = edown*edown/ydown;
    const Double_t ndown = ydown/scaledown;

    const Double_t eff = nup/ndown;
    if(eff>1){
      printf("ProperDivide ndown<nup!! %e %e %e\n", ndown, nup, eff);exit(1);
    }
    
    const Double_t sigma = sqrt(eff*(1-eff)/ndown);

    const Double_t rr = yup/ydown;
    const Double_t er = sigma*scaleup/scaledown;
    hup->SetBinContent(ii, rr);
    hup->SetBinError(ii,er);
  }

}

void drawMCbins0(const TString fin)
{
  TFile *ff=new TFile(fin);
  if(!ff){
    printf("no fin! %s\n", fin.Data()); exit(1);
  }

  TList *ll=new TList;

  TH2D * hh = (TH2D*)ff->Get("kFailEventVsPt");
  if(!hh){
    printf("no hh!\n"); exit(1);
  }
  //rebin here!!
  //has to be 4, otherwise if not an exact dividor, the the overflow will come and blow up the last bin.
  hh->RebinX(4);

  ll->Add(hh);

  //===

  TH1D * Psh=hh->ProjectionX("Psh",0,10);
  ll->Add(Psh);

  TH1D * FailPrimh=hh->ProjectionX("FailPrimh",2,2);
  if(FailPrimh->Integral(0,FailPrimh->GetNbinsX()+1)>1e-6){
    printf("bad FailPrimh! %f\n", FailPrimh->Integral(0,FailPrimh->GetNbinsX()+1));exit(1);
  }

  TH1D * FailVrech=hh->ProjectionX("FailVrech",3,3);
  ll->Add(FailVrech);

  TH1D * FailVzh=hh->ProjectionX("FailVzh",4,4);
  ll->Add(FailVzh);

  TH1D * PsVrecVzh=hh->ProjectionX("PsVrecVzh",1,1);
  ll->Add(PsVrecVzh);

  TH1D * PsVrech = (TH1D*)PsVrecVzh->Clone("PsVrech");
  PsVrech->Add(FailVzh);
  ll->Add(PsVrech);

  //==
  normal(Psh, PsN);
  normal(FailPrimh);
  normal(FailVrech);
  normal(FailVzh);
  normal(PsVrecVzh, PsVrecVzN);
  normal(PsVrech, PsVrecN);

  //===
  TH1D * hbin0Vz = (TH1D*)PsVrecVzh->Clone("bin0Vz");
  ll->Add(hbin0Vz);
  ProperDivide(hbin0Vz, PsVrech);
  hbin0Vz->Fit("pol0");
  const Double_t rbin0Vz=hbin0Vz->GetFunction("pol0")->GetParameter(0);
  const Double_t ratioVz = PsVrecN/PsVrecVzN;
  const Double_t errRatioVz = ratioVz/rbin0Vz - 1;
  printf("bin0resultA %s rbin0Vz  %e ratioVz %e errRatioVz %e\n", fTag.Data(), rbin0Vz, ratioVz, errRatioVz);

  TH1D *hbin0Vrec = (TH1D*)PsVrech->Clone("bin0Vrec");
  ll->Add(hbin0Vrec);
  ProperDivide(hbin0Vrec, Psh);
  hbin0Vrec->Fit("pol0");
  const Double_t rbin0Vrec=hbin0Vrec->GetFunction("pol0")->GetParameter(0);
  const Double_t ratioVrec = PsN/PsVrecN;
  const Double_t errRatioVrec = ratioVrec/rbin0Vrec - 1;
  printf("bin0resultB %s rbin0Vrec %e ratioVrec %e errRatioVrec %e\n", fTag.Data(), rbin0Vrec, ratioVrec, errRatioVrec);

  //===
  const Double_t PsVrecVzCorrN = PsVrecVzN * rbin0Vz * rbin0Vrec;
  TH1D * PsVrecVzCorrh=hh->ProjectionX("PsVrecVzCorrh",1,1);
  ll->Add(PsVrecVzCorrh);
  normal(PsVrecVzCorrh, PsVrecVzCorrN);

  //===

  TFile *fout=new TFile(Form("bin0%s.root", fTag.Data()),"recreate");
  ll->Write();
  fout->Save();
  fout->Close();

  printf("done!\n");
}

int main(int argc, char * argv[])
{
  for(Int_t ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  if(argc!=2){
    printf("argc!=2"); return 1;
  }

  ini(argv[1]);

  printf("fTag %s, PsN %.0f PsVrecN %.0f PsVrecVzN %.0f\n", fTag.Data(), PsN, PsVrecN, PsVrecVzN);

  drawMCbins0(argv[1]);

  return 0;
}
